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THE CALCULATION OF THE EIGENVALUES AND 
EIGENFUNCTIONS OF MATHIEU'S EQUATION 


I. INTRODUCTION 

Mathieu's equation, Eq. (1)> arises when the scalar Helmholtz 
equation is solved in the elliptic cylinder coordinate system by 
means of separation of variables. 

2 

(1) (1 - z 2 ) _ z g+ ( a + 2q - 4qz 2 ) f = 0 


Here the independent variable, z, may be related to either the angular 
or the radial elliptic variable, q is related to the ell ipti city of 
the coordinate system, and a is the separation constant or eigenvalue 
(also often refered to as the characteristic value). The solutions 
of this equation have been discussed extensively in the literature, 
and a summary of this literature may be found in Reference 1. The 
elliptic coordinate system is also discussed in detail in Reference 2. 

In the following a simple, direct method for the calculation 
of the eigenvalues and eigenfunctions of this equation is developed. 
The computational procedure is both rapid and accurate. This method 
is quite similar to a technique presented earlier[3] for the calcu- 
lation of the eigenvalues and eigenfunctions of the oblate and prolate 
spheroidal wave equations. 


II. THE EIGENVALUES 

The r 1 " solution of Eq. (1) may be expanded in terms of the 
trigonometric functions. 


( 2 ) 


f r (q»z) = 


, y 

fZ m=l 


[Ajj(q) cos (mz) + B^(q) sin (mz)] 


Note that the normalization of the leading expansion coefficient differs 
slightly from that customarily found in the literature. If this ex- 
pansion is substituted into Eg. (1) four independent sets of recursion 
relations are obtained, Eqs. (3-1 1 ) 


1 



Case 1 : 

(3) aA Q -y?qA 2 = 0 

(4) - /2- q A q + (a-2 2 ) A g - q A 4 = 0 

(5) - q A m _ 2 + (a-m 2 ) A m - q A m+2 = 0, m = 4, 6, 8, 


Case 2: 

(6) (a-l 2 -q) A 1 - q A 3 = 0 

(7) - q A m _ 2 + (a-m ) A^ - q A m+2 = 0, m = 3, 5, 7,•• , 


Case 3: 

(8) (a-2 2 ) B 2 - q B 4 = 0 

(9) " q B m-2 + ( a ' m2 ) B m " q B m+2 = °> m = 4, 6, 8, ••• 

Case 4: 

(10) (a-l 2 +q) B ] - q B 3 = 0 

01) " q B m _ 2 + (a-m 2 ) B m - q B^ = 0, m = 3, 5, 7, ••• 

Here both the eigenvalue, a, and the expansion coefficients, A m and 
B m , are dependent upon the order, r, of the solution; however, this 
dependency has been temporarily suppressed in the notation for the 
sake of simplicity. The choice of the particular set of recursion 
relations to be used is dependent upon the symmetry and periodicity 
of the desired solution. It is evident that the recursion relations 
for even solutions will involve only the A£ coefficients while the 
recursion relations for odd solutions will involve only the B£ coef- 
ficients. Further, if the solution has tt periodicity in z then only 
coefficients having even subscripts, m, will be utilized. Similarly, 
if the solution has 2tt periodicity in z then only coefficients having 
odd subscripts, m, will be utilized. These conditions are summarized 
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in Table I along with a shorthand notation which will be used to denote 
these cases. The integers s and t will be used to identify the 
periodicity and symmetry of the desired solution as indicated in Table I. 
In addition the letter D is used to denote the appropriate A or B 
coefficient as shown in Table I. The subscripts have been shifted so 
that D 0 is the leading coefficient in each case; this permits some 
simplification in subsequent programming. 


■ 

PERIODICITY 

TT 

s = 0 

2ir 

s = 1 

SYMMETRY 

EVEN 

t=0 

a£ j* 0, m even 
m 

^ = A m 

m m 

a![ f 0, m odd 
m 

D r = A r 

m Vl 

O i— 

o II 
O -P 

b][ f 0, m even 
m 

D m = 

f 0, m odd 

n r - R r 
m B m+1 


TABLE I 

SYMMETRY AND PERIODICITY CONDITIONS 


The adoption of the notation described above allows the four sets 
of recursion relations to be written compactly as one statement. 


(12) - + <VV q » D m (q) ' V C2 (q) * °« 

for m = 0, 2, 4,* •• 


where 


(13) X m = /2 if s = 0, t = 0, m = 0 

= 1 otherwise, 

(H) W m = [m + s + 2 t(l-s)] 2 + V m , 
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and 


V m = + q ifs = l,t = 0, m = 0 

(15) = - q if s = 1, t = 1, m = 0 

= 0 otherwise. 

The set of equations defined by Eq. (12) may be written in matrix form 
as shown in Eq. (16) 



The resulting real matrix is both symmetrical and tridiagonal. If 
this matrix is truncated to an N-by-N matrix, the bisection method[4] 
may be used to determine the eigenvalues, a~, in a rapid, accurate manner. 
This procedure has the distinct advantage tnat the speed of convergence 
is known prior to computation, i.e., the uncertainty of the unknown 
eigenvalue decreases by a factor of 2 upon each iteration. Test 
calculations have indicated that the truncation of the matrix to a 
dimension only slightly larger than the order of the largest eigenvalue 
produces no significant error. The N eigenvalues determined in this 
manner are denoted by a r , r = s, s + 2, s + 4, ♦••,$+ 2N-2, when 
arranged in order of increasing algebraic value. The bisection pro- 
cedure is initiated by noting that 
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This procedure has been tested for various cases with 0.05 _< q ^5.0 
and r max up to 30. In all test cases the bisection procedure was 
terminated after obtaining five significant figures. A sample set of 
eigenvalues is given in Table II. The computed eigenvalues were 
found to agree with tabulated values. [1,5] The listing of the sub- 
routine for the eigenvalues is given in Appendix A. 


ORDER 

Hodge 

Abramowitz & Stegun 

0 

-0.5800066E+01 

-0.580004602E+01 

2 

0. 74491 00E+01 

0. 74491 0974E+01 

4 

0.1709660E+02 


6 

0.3636098E+02 


8 

0.6419890E+02 


10 

0.1001261E+03 

0 . 1 001 2636922E+03 

12 

0.1440878E+03 


14 

0.1960645E+03 


16 

0.2560488E+03 


18 

0.3240378E+03 


20 

0.4000327E+03 


22 

0.4840264E+03 


24 

0.5760225E+03 


26 

0.6760208E+03 


28 

0. 78401 46E+03 



TABLE II 

Sample set of eigenvalues for q = 5.0, even symmetry, 
and it periodicity. (5 significant figures) 


III. THE EIGENFUNCTION EXPANSION COEFFICIENTS 

Having obtained the eigenvalue, a r , the eigenfunction expansion 
coefficient, D^, are obtained readily by means of recursion. Since 
these coefficients reach their maximum value at r % m, the recursion is 
carried out in two directions. Equation (12) is used to recur upward 
until m = r. 


07 ) 



< a r - V 

x m q 
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where 


D*" = 1 and D r , = 0. 
o -2 


Similarly, Eq. (12) is used to recur downward from m = rn v to m = r, 

ITlaX 


(18) 







where D m = 10”^°, 
max 


Dm max+2 



and m max is taken to be significantly 


larger than r. Values of as small as r + 5 have been used 
successfully; however, this choice will depend upon q and the accuracy 
desired. The set of coefficients for r > m are then normalized such 
that the two sets agree for m = r in order to obtain one consistent set 
of expansion coefficients. 


Finally the normalization condition 


(19) 


7T 

f 2 (z,q) dz = ir/2 
o 


is imposed. Substituting Eq. (2) into Eq. (19) yields 


( 20 ) 




Eigenfunction expansion coefficients obtained in this manner have 
been computed for the same ranges as the eigenvalues discussed earlier 
with i% ax up to 40. The agreement with tabulated values in all cases 
where possible was excellent. The listing of the subroutine for the 
eigenfunction expansion coefficient is given in Appendix B. 

The program utilizing these procedures required a total compilation 
and execution time of 12 seconds on an IBM 360-75 computer for the 
calculation of 60 eigenvalues and 1,260 eigenfunction expansion coef- 
ficients. 
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IV. CONCLUSION 

The eigenvalue problem associated with Mathieu's equation when 
cast in matrix form was found to yield a real, symmetric, tridiagonal 
matrix. Thus the highly efficient and accurate bisection method may 
be used inmediately to determine the eigenvalues. The eigenfunction 
expansion coefficients are subsequently obtained by a standard recursive 
technique. This procedure was tested and the results were compared 
with previously tabulated results showing excellent agreement. No 
computational difficulties were encountered. 
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APPENDIX A 

EIGENVALUE SUBROUTINE 


The Fortran IV subroutine listed in Fig. A-l may be used to 
compute the first NMX eigenvalues, a (q), of Mathieu's equation. The 
input parameters are: 


NMAX = 2*NMX 

NMX = number of eigenvalues desired 
Q = ell ipti city 

IS = integer 0 or 1 depending upon the periodicity as defined 
in Table I of the text 

IT = integer 0 or 1 depending upon the symmetry as defined in 
Table I of the text. 


The NMX computed eigenvalues are output in the one-dimensional array, 
EIG. 


The number of significant figures in the computed eigenvalues 
is determined by the value of ACC which is defined near the beginning 
of the subroutine. In the case shown approximately five significant 
figures will be produced. If this accuracy is not obtained within 
60 iterations the subroutine will output the error statement shown in 
statement 42. The maximum number of allowed iterations may be adjusted 
by altering the constant in the statement preceeding statement 41. 
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SUbRGO I INE MA TF l G ( NMAX » 0 » I S » IT ,EIG) 

U I PENSION PIbO), IP(bO), ALPHA (t>0)» BETA (bO) ,ElG(t>0) 

ACOl.OL-05 

N=NNAX+2 

N I = N+ 1 

P( 1 ) = 1 .0 

1P< l) = l 

X=L 

1 F ( ( I S. EU. 0 ) • AND* ( IT. EG. 0)) X= 1.41421 36 
V=0 

I F ( CIS . EC. I ) • ANO. ( II .EQ.O) ) V=0 
I F ( (IS.EG.l). ANO. (IT.EQ.l)) V=-Q 

00 2 I — 1 -» N 

ALPHA { l )=- (2*( I-I)+IS+2*IT*( 1-IS) )**2-V 
bfcT A( l + l )=-X*« 

X=1 

2 V=G 

BETA(N+l )=0.0 

BO=AbS ( ALPHA ( l ) ) +ABS ! BETA I 2 ) ) 

DO 3 1=2, N 

A0= ABS l BET All) ) + ABS( ALPHA! I) ) +ABS( BETA! 1 + 1 ) ) 

BETA! I )=BETA( I )*BETA( I ) 

1 F ( AO- BO ) 3 , 3 * 4 
4 bO= AO 

3 CONTINUE 
AO=-BO 

600 NMX=NMA X/2 

DO 20 K= 1 » NHX 

A=AO 

B=BO 

IERR=-1 

21 IlS=0 
C0=(A+B)/2. 

IF(CO)50,22,60 

60 £RK= ( B-A ) /ABS ( CO ) 

I ERR= l ERR+l 
I F ( IERH-60)40,41,4l 

41 WRITE! 6,42 ) K 

42 FORMAT ( 3bH0 I ThRAT IONS EXCEEDED FOR EIGENVALUE ,13) 
GO TO 700 

40 l F ( ERR- ACC )24, 24,22 

22 P(2)=ALPHAU)-C0 
P( l ) = l .0 

00 5 1=3, NI 
ABP=ABS ( P ( 1-1 ) ) 

1 F ( ABP. L I . 10.0 ) Gfl TO b 
P(I-2)=P( I-2J/ABP 

P( l-l )=P( 1-D/ABP 


Fig. A-l . Eigenvalue subroutine. 
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5 P([ )=( ALPHA d-l> -CO) *P(I-1 )-BETA( I-l)*Pl 1-2) 
DO 6 I -2 » N1 

IF(PU) )1A»8*9 

8 IF(P(I-l))9,9,14 
l A IP(I)=-1 

GO TO 10 

9 1 P 1 1 1=1 

10 lF(IPd )-IP(l-l) )6,ll,6 

11 IIS=IIS+1 

6 CONTINUE 

IF( I1S-K) 16*15.15 

15 A=CO 

GO TO 21 

16 B-CO 

GO TO 21 
2 A BO=CO 
700 EIG(K)=-CO 
20 CONTINUE 
801 RETURN 
END 


Fig. A-l . (Contd.) 



APPENDIX B 

EXPANSION COEFFICIENT SUBROUTINE 


The Fortran IV subroutine for the eigenfunction expansion 
coefficients is listed in Fig. B-l . In addition to the parameters 
defined in Appendix A, this subroutine requires one additional input 
parameter, MAXR, which is the number of terms to be retained in the 
expansion. The expansion coefficients are output in the two- 
dimensional array, D(N,J), where N denotes the N™ eigenfunction and 
J denotes the Jtn expansion coefficient. 
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SUBROUTINE MA TCOF ( NM AX * 0 * I S • l T* MAXR, E IG* D ) 

DIMENSION EIG ( 50 ) *D ( SO* 50 ) *DP ( 50 ) 

NMX=NMAX/2 
DO l N=L»NMX 
DP(MAXR+3)=0. 

DP (MAXR+2 )=1 »0E— 30 
DIN, l)=0. 

D(N,2)=1.0 

X=1 

Y=1 

I F ( (IS.EQ.O). AND. ( IT. EQ.O I ) X=U4142136 
V=0 

IF (( IS.EQ.U.ANO.UT.EQ.OI) V=Q 
IF UIS.EO.D.ANO.IIT.EO.m V=-Q 
DO 107 LL=l *MAXR 
L=LL 

IF(LL.GT.N) GO TO 300 

D(N,L+2)=-(Y*D(N»m/X+IEIGCN)-I2*(L-l) + IS+2*IT*( 1-IS) ) 
l**2-V)*D(N,L+l)/(X*Q) 

Y=X 

X=1 

V=0 

GO TO 107 

300 L=MAXR+N-LL+l 

DP(L+l)=-DP(L+3)+(ElG(N>-(2*L+IS+2*IT*(l-IS> >**2-V) 
l*DP(L+2)/Q 
107 CONTINUE 

CONJOIN, N+2)/UP(N+2) 

DO 118 J=N* MAXR 
118 D ( N* J+2 )=CON*DP ( J+2 ) 

SUM-0 

MRX=MAXR+2 
DO 301 J=2 » MRX 

301 SUM=SUM+D ( N* J ) **2 
ALF-SQRT ( SUM) 

DO 302 J-2 * MRX 
D(N,J-1)=D(N* JI/ALF 

302 CONTINUE 
1 CONTINUE 

RETURN 

END 


Fig. B-l . Eigenfunction expansion coefficient subroutine. 
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